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We present local numerical models of accretion disk turbulence driven by the magnetorotational instability with varying 
shear rate. The resulting turbulent stresses are compared with predictions of a closure model in which triple correlations 
are modelled in terms of quadratic correlations. This local model uses five nondimensional parameters to describe the 
properties of the flow. We attempt to determine these closure parameters for our simulations and find that the model 
does produce qualitatively correct behaviour. In addition, we present results concerning the shear rate dependency of the 
magnetic to kinetic energy ratio. We find both the turbulent stress ratio and the total stress to be strongly dependent on the 
shear rate. 



1 Introduction 

Since the work of Balbus & Hawley (1991) it is now gen- 
erally accepted that turbulence in accretion disks is caused 
by the linear magneto-rotational instability (hereafter MRI) 
that was originally discovered by Velikhov (1959) in con- 
nection with Couette flow of liquid metals. This linear in- 
stability can be excited in sufficiently ionized, differentially 
rotating fluids where the angular momentum decreases with 
increasing radius. When the fluid is threaded by a mag- 
netic field, differential rotation causes stretching of the field 
lines. The tension that builds up opposes the shearing, act- 
ing to enforce rigid rotation. If the field is subthermal, the 
interplay between differential rotation and magnetic tension 
destabilizes the fluid, resulting in linear growth of Reynolds 
and Maxwell stresses, which transport angular momentum 
outwards. The instability eventually leads to a fully turbu- 
lent, non-linear state. 

Let us assume a rotation profile of the form ft oc r~'?, 
where is the angular velocity at distance r. A necessary 
condition for the MRI to be excited is (j > 0, i.e. the an- 
gular velocity decreases outward. The Keplerian case with 
(7 = 1.5 has been extensively studied by means of numerical 
simulations making use of the shearing box approximation 
(e.g., Hawley, Gammie & Balbus 1995, 1996; Brandenburg 
et al. 1995; Johansen & Klahr 2005), as well as in global 
disks (e.g. Armitage 1998; Hawley 2001; Lyra et al. 2008). 
These studies have shown that the MRI leads to a fully tur- 
bulent saturated state in which the Maxwell stress is respon- 
sible for the majority (about 80%) of the angular momentum 
transport. Even in the absence of vertical density stratifi- 
cation a turbulent small-scale magnetic field can be main- 
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tained by dynamo action (Hawley et al. 1996). When strati- 
fication is present, cyclic large-scale dynamo action can be 
excited (Brandenburg et al. 1995). 

The dependence of Reynolds and Maxwell stresses on 
the shear rate has been investigated numerically by Abramo- 
wicz, Brandenburg & Lasota (1996) in the presence of strat- 
ification and Hawley, Balbus & Winters (1999) in the ab- 
sence of stratification. It turns out that near q = the 
Reynolds stress, which couples to the large-scale vorticity 
W = {2 — q)il, is small due to the strong stabilizing ef- 
fect of the vorticity (Hawley et al. 1999). This makes the 
Maxwell to Reynolds stress ratio very high. As the shear 
rate q^l increases, the shear-coupled Maxwell and Reynolds 
stresses increase due to the decrease of the vorticity. How- 
ever, the growth of the Reynolds stress is significantly faster, 
so the stress ratio diminishes with increasing q. 

The interest in MRI-generated turbulent stresses as a 
function of q was rekindled in a recent study where a lin- 
ear analysis of the Reynolds and Maxwell stresses was pre- 
sented (Pessah, Chan & Psaltis 2006a, hereafter PCP06). 
They derived a simple relation for the stress ratio, depending 
only on the shear parameter q. Comparing this result with 
the non-linear simulations of Hawley et al. (1999), they find 
that, even in the saturated turbulent regime, the stress ra- 
tio does indeed depend almost entirely on q alone, and that 
there is only a weak dependence on other properties of the 
flow or on the initial conditions. 

The Shakura-Sunyaev viscosity parameter a (Shakura 
& Sunyaev 1973) is still a popular tool to link accretion 
disk observations to theory. The so-called a model of disks 
is based on the assumption that the turbulent stresses Tr^ 
scale linearly with the thermal pressure; a < 1 being the 
proportionality factor. In this formalism, the outward trans- 
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port of angular momentum is characterized by a turbulent 
viscosity of the form 

ft = aCsH, (1) 

where the eddy viscosity is assumed to scale with the sound 
speed Cg and the correlation length with the disk scale height 
H. 

As this model linearizes the turbulent system into a stan- 
dard Navier-Stokes fluid suitable for analytical manipula- 
tion, it has been an invaluable tool in developing theory 
of accretion disk dynamics, even though the parametriza- 
tion offers no explanation of what is causing the viscosity. 
Even after the realization of the astrophysical importance 
and subsequent reinvigorated interest in the theory of the 
MRI, the a-model is still widely used. 

New, more physically motivated models have been de- 
veloped recently to describe the angular momentum trans- 
port in accretion disks (e.g. Kato & Yoshizawa 1995; Pes- 
sah, Chan & Psaltis 2006b; Ogilvie 2003, hereafter O03). 
A common characteristic of these models is the treatment 
of the problem of angular momentum transport: the govern- 
ing magnetic and kinetic equations are divided into linear 
and non-Unear terms. The non-linear correlation functions 
are then modelled by closing the system of equations with 
approximate expressions that still embody the physics but 
are considerably easier to solve (such closure models are 
commonly used in modelling turbulence; see references in 
O03). It should be noted that the a-model is mathematically 
equivalent to the simplest closure model, where the correla- 
tion functions are modelled by one single coefficient. All 
these models also operate in the absence of mean magnetic 
fields. A first attempt to validate the O03 model was made 
by Garaud & Ogilvie (2005) in connection with shear flow 
where linear and non-linear instabilities were found to be 
well reproduced by the model. 

In the present study, the turbulent stresses are extracted 
as functions of q from non-stratified local numerical mod- 
els with zero net flux using the shearing box approxima- 
tion. The simulation data are compared with the linear re- 
sults of PCP06 and the non-linear closure model of O03. 
There are a few other closure models that describe the be- 
haviour of local, magnetohydrodynamic turbulence. In par- 
ticular, the model of Pessah et al. (2006b) assumes a uni- 
form vertical field in the disk which is not incorporated in 
our three-dimensional simulations. Therefore we concen- 
trate here mainly on the closure model of O03. 

The remainder of the paper is organised as follows: in 
Sections |2] |3] and |4] the linear MRI model of Pessah et al. 
(2006b), the non-linear closure model of O03, and the nu- 
merical model are presented, respectively. Furthermore, the 
results and related discussion are presented in Sections |5] 
and|6] 

2 Stress ratio in the model of PCP06 

Here we briefly summarize the formalism employed in the 
model of PCP06. Using the kinematic hydromagnetic equa- 



tions, PCP06 calculated the ratio of the relevant components 
of Reynolds and Maxwell stresses in a local approximation. 
We use here a Cartesian frame of reference, where x, y and 
z denote the radial, toroidal and vertical directions, respec- 
tively. The relevant component of the stress tensor is then 
Txij, which can be decomposed as 



T 

where 



Mx 



(2) 



(3) 



are the Reynolds and Maxwell stresses, respectively, and u 
is the departure from the mean flow, b is the departure from 
the mean magnetic field, p is the density, po the vacuum 
permeability, and angular brackets denote a suitable volume 
average. Throughout our paper, compressibility effects are 
ignored in the turbulence models, i.e. p ^ po is assumed 
constant, even through the simulations are fully compress- 
ible. The rms velocity in our simulations remains consider- 
ably smaller than unity and the flow is therefore subsonic. 

Using linear theory, PCP06 found that for given wave- 
number k the stress ratio is given by 



Mxy{k) 



R: 



xy{k) 



= 1 



2(2-g)f]2 
Pi 



(4) 



where 57o is the angular velocity, wa = Bq/^pqPq is the 
Alfven speed based on a constant vertical magnetic field, 
and 7fe is the corresponding growth rate of the mode with 
wavenumber k. Moreover, for the fastest growing mode 
with 



7femax/f^0 

they find 



4-g 



2 7 2 /r>2 1 2 



Rx 



(5) 



(6) 



This expression shows that, for the relevant case with q < 2, 
the magnitude of the Maxwell stress is always larger than 
that of the Reynolds stress. This formula provides a strik- 
ingly simple prediction of the stress ratio which is in good 
agreement with simulation data, even if there is no imposed 
magnetic field (PCP06). One envisages that the relevant 
wavenumber is able to adjust itself to the value where the 
growth rate is maximal. 

3 The O03 closure model 

The local closure model of O03 includes the linear inter- 
action of the turbulent stress tensors with shear and rota- 
tion. The linearized evolution equations for the Maxwell 
and Reynolds stresses can be derived directly from the basic 
MHD equations and are fairly straightforward to solve nu- 
merically. For modelling the non-linear triple correlations 
of fluctuating quantities and the small-scale diffusion, five 
dimensionless, positive definite coefficients appear in the 
closed system of equations. A closure model is needed to 
deal with these non-linear terms which are described by 
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physical effects and constrained by symmetry properties 
and dimensional considerations. 

The five closure coefficients stand for the turbulent dis- 
sipation of the Reynolds stresses (Ci), their isotropization 
(C2), the effect of the small-scale Lorentz-force as a source 
for Rij combined with a sink for (C3), a source due 
to small-scale dynamo for combined with a sink for 
Rij (C4), and turbulent dissipation of the Maxwell stresses 
(C5). 

The evolution equations of this model are given by 



dtRij 



-1^-1/2 



c. 



M 



-1 -1/2. rM 



(7) 

(8) 



where t = L/U is the turnover time, U = i?^/^ is the rms 
velocity, L is the typical scale of the energy-carrying eddies, 
and 



R = R^^^ {p{l 



M = M„ 



■bD/po, 



(9) 



(10) 



are the traces of tensors Rij and Af^ . Furthermore, C^j and 
Affj are linear and nonlinear terms, respectively, for a ~ R 
or Af, and are given by 

C^j = RikU j^k + RjkUi,k + '2^k{£jklRil+£iklRjl)lii) 

C^^ M,kUj,k + Mjk^^,k, (12) 
A/^^= {Ci + CiB^) % - CaSAfy , (13) 
C^B^R,, - (C3 + C5)BM,^ , (14) 

where B ~ (M/R)^/^ is the ratio of rms magnetic and 
velocity fields. 



C2{Rv 



\R5.,,) 



is an isotropization term, and Ci . . . C5 are positive con- 
stants that are of the order of unity. Advection operators of 
the form U kdk have been neglected. The contribution of 
the advection terms vanish under fully periodic boundary 
conditions on average so they make no contribution. Note 
also that t and B are time-dependent, because R and M are 
time-dependent. 

A similar model for the hydrodynamic case has recently 
been used to to model the generation of shear in rotating 
anisotropic turbulence by the A effect; see Kapyla & Bran- 
denburg (2008) for details. An important difference is that 
in their model the flow was driven by an external body force, 
so the equations for the evolution of the Reynolds stress 
have a corresponding forcing term as well. Such a term is 
here absent, because the turbulence is solely the result of 
shear flow instabilities and are modelled by the equations 
of 003 without external forcing. As explained by O03, the 
model also predicts turbulence for q <i), where simulations 
have not shown self-excited turbulence. The reason for this 
is that 003 do not specifically model the MRJ dynamics. 

In the present study the large scale velocity is given by 



the shear flow U ^ Un 



-qfloxey. We abbreviate the 



parameter combination Ci/L of O03 by c^. O03 gave an 



analytic solution for the hydrodynamic case. Here we give 
a partial solution for the components of the two stresses if 
R and M are known. In the steady state, the equations yield 
for the Reynolds stresses 



i?,, = RA (ic2 + Q) , 

Ryy = RA (ci + \c2 + aB^cs - Q) , 

Rzz = ^R-^ C2, 



R 



2qn 



[ciR - {Bc3 - C4)M] 



where we have introduced the abbreviations 



A = (ci + C2 + aS^cs) 



C4 



C3 + C5 ' 



and 



Q^-[ci-B^ {Bc3 - C4) 



(16) 
(17) 
(18) 

(19) 
(20) 
(21) 



The components of the Maxwell stress tensor can be ex- 
presses in terms of the corresponding components of the 
Reynolds stress tensor as 



'-xx — 0,BRxx 5 



Myy = aB {Ryy ^ R)+ M, 



Af. 



2qn 



[C4 - (C3 + C5)B] M. 



(22) 
(23) 
(24) 

(25) 



(15) The remaining mixed components involving z vanish, i.e. 



Mr 



A/„ 



0. When A/ = 0, there is an 



explicit expression for R, derived by O03. 

For the general case with A/ 7^ 0, we did not find an 
analytic expressions for R and A/ in closed form. However, 
using the time-dependent equations it is possible to deter- 
mine a linear fit of the form 



R = R^°'^ + J2 n U - a 



JO) 



A/ = A/("' + > TO, 



5 

1=1 



(26) 



(27) 



where the Cj"^g are an approximation to the final fit parame- 
ters, and i?(0) and A/("' are the corresponding numerically 
determined values of R and M. We are thus able to decrease 
the degree of freedom of the fitted system. We refine the pa- 
rameters Ci with respect to the initial "guess" cf^ based on 
the quantity 



5 = [(A^O03 - M,,^f + (i?O03 - i?shn)'] ^'^ ■ (28) 

In the analysis we seek the minimum value of &. 
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4 Simulations 

For the simulations we adopt a cubic computational domain 
of size (27r)''. The gas is isothermal with constant sound 
speed Cs', vertical gravity and stratification are omitted. The 
calculations are local and therefore made under the shearing 
box approximation (see below). The equations to be solved 
read 

^ = -V • {pu) , (29) 
Du 

— = -(tt • V)tt - qV,oUxey - 2 QqEz x u 

--Vp+-JxB + -V-(2iypS), (30) 

PR P 

VA 

— = ux B + qfloAye^ + rjV^A , (31) 

where V/Vt = d/dt + Uod/dy includes the advection 
by the shear flow, u is the departure from the mean flow 
Uq, p is the density, A is the magnetic vector potential, 
B = V X A is the magnetic field, and J ~ 'V x B/fiQ 
is the current density, po is the vacuum permeability, = 
^{ui,j + Uj,i) — ^Sij'V ■ u is the traceless rate of strain 
tensor, ly is the kinematic viscosity, and 77 is the magnetic 
diffusivity. 

In order to get as close as possible to the ideal limit, 
we replace the diffusion terms by a hyperviscosity scheme, 
i.e. we replace the V^-operators by V^, aiming at maxi- 
mizing the Reynolds number in the quiescent regions of the 
flow while diffusing and damping fluctuations near the grid 
scale. Compared to direct simulations (e.g. Haugen & Bran- 
denburg 2006) with uniform viscosities, smaller grid reso- 
lution can be used to resolve the flow, which is crucial if the 
plan is to, e.g., undertake a parameter study, like in our case. 

Periodic boundary conditions are applied in all three di- 
rections; in the radial direction we account for the shear flow 
Uq by making use of the shearing box approximation (e.g. 
Wisdom & Tremaine 1988): 

f{\L.j,,y, z) = f{-jL.j,,y + q^oL^t, z), (32) 

where / stands for any of the seven independent variables, 
Lx stands for the radial extent of the computational domain, 
and t is the time. 

The domain is initially threaded by a weak magnetic 
field, 

A ^ AoBy cos kx cos ky cos kz , (33) 

where k = ki has been chosen, and fci = 2tt/Lz is the 
smallest finite vertical wavenumber in the domain of height 
Lz. Thus, the magnetic field contains periodic x- and z- 
components with amplitude ylo- 

We choose the values of /ca, 51o and Aq so that the most 
unstable mode of the MRI, fcmax = ^ /ua is well resolved 
by the grid; in practise this means that we always adopt 
k/ki ^ 1, Qq = 0.2csfci and Aq ~ 0.2^pqPq Cgfc^^, re- 
sulting in fcniax — 0{ki). For the initial setup, the other 
condition for the onset of MRI, namely P ^ 1, where /? is 



1.0 




k v,/n 

Fig. 1 Growth rates as a function of wavenumber from the 
ID calculations. Crosses: q ~ 1.75, triangles: q — 1.5, dia- 
monds: q = 1.25, squares q = 1.0. The solid lines represent 
the linear growth rates; see, e.g., PCP06. 

the ratio of the thermal to magnetic pressure, is also satisfied 
as P is minimally 50 at the maximum values of the magnetic 
field. In all our runs the Mach number is well below unity, 
so compressibility effects are negligible. 

For all the simulations we use the Pencil-Cod^ 
which is a high-order (sixth order in space, third order in 
time), finite-difference code for solving the MHD equations 
(Brandenburg & Dobler 2002). 

The local calculations have been carried out at two dif- 
ferent resolutions, namely 128 (in ID) and 64'^ (in 3D); the 
corresponding numerical diffusion coefficients are z^hypcr = 
'7hypcr = 3.5 X 10^^ and 2.0 x 10^^, respectively. The cal- 
culations were carried out on the IBM eServer Cluster 1600 
at Scientific Computing Ltd., Espoo, Finland. 

5 Results 

5.1 Linear results 

In order to make contact with the model of PCP06, we first 
aim at reproducing their theoretical linear results (see also 
Balbus & Hawley 1991) using one-dimensional calculations 
with an imposed magnetic field. We calculate several sets by 
fixing the shear rate {q = 1.00, 1.25, 1.50, 1.75) and angular 
velocity ilo and varying the initial magnetic field strength 
Bq. The growth rate and the stress ratio are monitored dur- 
ing the exponential growth of the instability. The results are 
displayed in Figs.[T]and|2l 

As can be seen from Fig. [T] the linear growth rates can 
be reproduced by the numerical method quite accurately. 
From Fig. |2] it can be observed that the linear prediction 
(dotted curve) intercepts the numerical results (solid lines) 
exactly at fcmax for each q-curve. At each fcmax (<z), there- 
fore, the linear prediction and numerical results show per- 
fect agreement; see Eq. In the ID calculations, however, 

' http://www.nordita.org/software/pencil-code/ 
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Fig. 2 Stress ratio ~Mxy/Rxy as function of wavenum- 
ber. The solid lines represent the linear ID results with vary- 
ing shear parameter q. The dotted line shows the linear pre- 
diction of PCP06 plotted for each fcmax('?)- 

the system has not much freedom to create any other MRI 
mode than the one that is imprinted by the initial magnetic 
field strength, due to which the stress ratio is observed to 
vary as function of wavenumber, so that monotonically de- 
creasing stress ratios are found with increasing wavenum- 
ber for all the qs investigated. This is clearly in disagree- 
ment with the PCP06 assumption, according to which the 
mode with fcmax should always get preferentially excited. 
The wavenumber dependence of the stress ratios, however, 
disappears in 3D: independent of the initial magnetic field 
strength, the mode with fcmax is observed to dominate. In 
that sense our 3D results are giving support to the basic as- 
sumption of PCP06, although, as will be discussed in the re- 
maining part of the paper, the magnitudes and (/-dependence 
of the stresses are otherwise different from the linear analy- 
sis. 



5.2 Nonlinear results 

We have performed a set of local calculations, in which 
we have fixed the angular velocity JIq and the strength of 
the initial magnetic field, and then varied the shear pa- 
rameter q. For the investigated range of shear parameters 
0.4 < q < 1.9, the maximally growing wavenumber, there- 
fore, is varying according to Eq. (|5]), but is resolved by the 
numerical grid in all the calculations. We have checked that 
the resulting stresses are independent of the initial magnetic 
field strength (computations with initial plasma (3 of 20- 
800 were performed). This is particularly important in the 
net flux case (Blackman, Penna & Varniere 2008). 

For smaller shear the stabilizing effect of vorticity is 
strong and the fluctuations grow larger than the stresses 
themselves; these calculations, therefore, are not included 
in the results. The data averaged over the nonlinear satu- 
rated stage (typically from a hundred to a few hundreds of 



rotations) for the range 0.4 < q < 1.9 is presented in Ta- 
ble □ 

As is evident from Table [T] both Reynolds and Maxwell 
stresses grow with the shear parameter q. The growth of 
the Reynolds stress is much stronger than the growth of the 
Maxwell stress, due to which the stress ratio (plotted in the 
second panel of Fig.|3]l decreases as function of q. The simu- 
lated stress ratio significantly differs from the PCP06 linear 
prediction given by Eq. Q and plotted in the figure with 
dotted lines. All the data points including the error bars are 
consistently larger than the prediction by a factor of 2-3. 

Panel 1 of Fig. |3] shows the total stress, defined in 
Eq. (|2]i, as a function of the shear-to-vorticity ratio, q/ (2 — 
q). The shear-to-vorticity ratio is a quantity that is inde- 
pendent of the coordinate system, which was the main rea- 
son why Abramowicz et al. (1996) presented the stress as a 
function of this ratio. Interestingly, they found that the stress 
is a nearly linear function of the shear-to-vorticity ratio. This 
is confirmed by the new data. If one were to plot the stress as 
a function of q directly, the relation would become strongly 
nonlinear 

Also the magnetic to kinetic energy ratio (panel 6 in 
Fig. |3]l exhibits q dependency, but it is less strong than that 
of the magnetic to kinetic stress ratio (panel 2 in Fig.O. The 
energies are defined as Ek — ^(pu^) and Em = ^(B^). 
Thus Ek = ii? and Em = 5^^- Our results indicate that 
for flat (galactic) rotation curves with q = 1 the energy ratio 
should be a factor of two higher than for the case of Keple- 
rian accretion disks (q = 1.5). 

5.3 Predicting stresses with the O03 model 

In their treatment, O03 use fiducial closure parameters 
Ci_,5 = 1 in order to demonstrate the overall behavior of 
the model. He recommends the parameters to be calibrated 
by comparison with numerical simulations in order to ob- 
tain more accurate predictions. We have made an attempt to 
determine the dimensionless closure parameters that work 
for each value of q in the range 0.4 < q < 1.9. 

There are essentially two ways to approach the problem. 
The first is to take the time independent equations (fT6l)-(l25l) 
and solve for ci_5 using values of the traces R and M from 
the simulations. We call this method "backward modelling". 
The problem with this method is that the Eqs. (fT6])-(l25]) are 
incomplete and do not set any constraints to the ratio AI / R. 
Consequently, it does not produce the correct solution to 
the time dependent equations (|7]i and ([S]). Therefore we also 
used the linear approximation described by Eqs. ( l26b and 
( l27b to improve the fit. 

The second way of determining ci_5 is what we call 
"forward modelling". The idea here is to seek such ci_5 
that the results of the time-dependent equations (|7| and ([8]) 
yield the same individual stresses and the traces R and M 
as the numerical simulations. This time a universal ci_5 
is determined so that it predicts the stresses for all values 
0.4 < q < 1.9. Once a reasonably good set of ci_5 was 
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Table 1 Stress components averaged over the saturated regime of the local calculations. For the O03 closure model, 

R'xz = Ryz = = Myz = 0. 
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Fig. 4 The error estimate S (see Eq. (28)) as a function of 
C. 

found, the result was fine tuned further using the linear ap- 
proximation given by Eqs. ( |26] l and ( |27l ); see Figs.|4]and|5] 
From the Fig.|4]it can be seen that the best fit is obtained by 
fixing C4 while keeping the other Ci unchanged. 

The final closure parameters are thus ci = 0.63, C2 = 
0.73, C3 ~ 0.33, C4 = 0.58 and C5 — 1.35, corresponding 
to d = c,L with C\ ^ 4.0, C2 4.6, C3 = 2.1, C4 = 3.6 
and C5 = 8.5. For these ci_5 the relevant fit parameters in 
Eqs. (EUl and ^ are = 0.0057 and il/t") 0.013 
together with ri = -0.0065, r2 = -0.00012, = 0.024, 
r4 = -0.010, rs = -0.0061, and mi -0.0057, mz = 
-0.0012, ma = 0.0042, TO4 = 0.0038, mg = -0.018. 

6 Discussion 

In this study we set out to investigate the shear rate depen- 
dency of MRI-generated turbulent stresses. We have per- 
formed a series of local shearing box simulations with vary- 
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Fig. 5 Comparison of the simulation data (diamonds) to 
the chosen fit to Ogilvie model with ci = 0.63, C2 = 0.73, 
C3 = 0.33, C4 = 0.58, C5 = 1.36 (diamonds) for one 
particular run with q=l.5. The solid line shows the M/R- 
dependent result of the stationary equations with the chosen 
set of c-parameters. The additional constraint used in the fit- 
ting procedure required that the difference in the simulated 
and model AI/ R is minimal, due to which in this solution 
the values of M/R match while the individual stress ratios 
somewhat differ. 



ing q and measured the resulting turbulent stresses. We 
find that the turbulent stress ratio —AIxy/ Rxy, and the to- 
tal stress Txy exhibit strong (/-dependency. The relation for 
the stress ratio by PCP06; see Eq. Q predicts similar be- 
haviour, but the ratio computed from the simulations is con- 
sistently 2-3 times larger than what their result indicates. 

In order to further study the evolution of the MRI and 
the stresses we have attempted to reproduce them using the 
local closure model by O03. We first find a set ci_5 such 
that the time-dependent equations give the same individ- 
ual stresses and traces R and M. The linear approxima- 
tion of Eqs. (I26b-(l27b is then used together with the time- 
independent equations to improve the fit. The O03 closure 
parameters that describe our numerical simulation results 
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Fig. 3 Panel 1: total stress as a function of the shear-to-vorticity ratio, q/ (2 — q), from the simulations (diamonds) and 
prediction from the O03 closure model (line) using ci = 0.63, C2 = 0.73, C3 = 0.33, C4 = 0.58 and C5 ~ 1.35. The 
error bars show the standard deviation of the quantity in question. Behaviour predicted by the Shakura-Sunyaev a viscosity 
model has been presented for comparison. Panel 2: stress ratio as a function of q from the simulations (diamonds), over- 
plotted with the 003 closure result (solid line) and the PCP06 linear prediction (dotted line). Panels 3-5 show several stress 
component ratios from the simulations overplotted with the O03 closure results. To help visualisation, stress components 
Rxx, Mxx have been scaled up with a factor of 6. and Ryy, Myy with a factor of 20. In panel 6, the ratio of magnetic to 
kinetic energy is presented. Throughout the figure, (s) is used to denote a simulation result and (c) a result given by the 
closure model. 



are found to be Ci = 4.0, C2 = 4.6, C3 = 2.1, C4 = 3.6 
and C5 = 8.5. 

The closure model by O03 was thus found to predict 
our simulation generated turbulent stresses quite well. The 
model certainly offers a much-needed method for studying 
the evolution of turbulent stresses in the shearing sheet limit. 
However, additional constraints for R and M are needed 
in order to determine the closure parameters Q. The linear 
fitting approximation described in Eqs. ( |26] ) and ( |27] | was 
devised to overcome this shortcoming. 

Our results may also have some relevance in the galactic 
context, in which the MRI has been proposed to be respon- 
sible of the anomalous turbulent velocity dispersions found 
in the outer regions of some galaxies (e.g. Sellwood & Bal- 
bus 1999). On the other hand, the energy balance estimates 
from observations of NGC6949 (Beck 2004) indicate that 
magnetic energy could become dominant over the kinetic 



energy in the outer regions of this galaxy, so that the en- 
ergy ratio Em j Ek ~ 3-4. According to our present results, 
in MRI-driven systems the magnetic energy clearly domi- 
nates over the kinetic energy for all g < 1.6, at the galactic 
value of g = 1 for flat rotation curves we find the value 
w 3, a number agreeing rather well with the observational 
estimates. 
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